Problème torsion d’arbres#
Librairies#
Show code cell content
# Importation des librairies
import dolfinx # Module principal de FEniCSx pour le calcul par éléments finis
from dolfinx import mesh, fem, plot, io, default_scalar_type # Sous-modules FEniCSx essentiels
from dolfinx.fem.petsc import LinearProblem # Résolution de systèmes linéaires
from mpi4py import MPI # Interface pour le calcul parallèle
import ufl # Unified Form Language pour les formulations variationnelles
import numpy as np # Calcul numérique efficace
import matplotlib.pyplot as plt
import math # Module de fonctions mathématiques standard de Python
import pyvista as pv # Visualisation 3D scientifique
import gmsh # Générateur de maillage 3D
import meshio # Lecture/écriture de différents formats de maillage
from dolfinx.io import XDMFFile # Gestion des fichiers XDMF pour données volumineuses
import dolfinx.fem as fem # Module FEniCSx pour la méthode des éléments finis
import dolfinx.mesh as mesh # Module FEniCSx pour la gestion avancée des maillages
from petsc4py import PETSc # Suite de solveurs numériques pour problèmes à grande échelle
import panel as pn
from myst_nb import glue
Géométrie#
Définition de la géométrie et du maillage#
Show code cell content
# Initialisation de GMSH
gmsh.initialize()
gmsh.model.add("ellipse")
# Définition des paramètres géométriques et de maillage
a = 0.2 # Demi-grand axe de l'ellipse en mètres
b = 0.2 # Demi-petit axe de l'ellipse en mètres
h = 1.0 # Hauteur du cylindre elliptique en mètres
# Calcul du rayon équivalent du cylindre elliptique
R = 1/2 * np.sqrt(a**2 + b**2)
# Création de la géométrie : ellipse de base
ellipse = gmsh.model.occ.addEllipse(0, 0, 0, a, b)
# Création d'une boucle fermée à partir de l'ellipse
curve_loop = gmsh.model.occ.addCurveLoop([ellipse])
# Création d'une surface plane à partir de la boucle
surface = gmsh.model.occ.addPlaneSurface([curve_loop])
# Extrusion de la surface pour créer un cylindre elliptique
gmsh.model.occ.extrude([(2, surface)], 0, 0, h)
# Synchronisation du modèle géométrique
gmsh.model.occ.synchronize()
# Configuration des paramètres de maillage
lc = 0.02 # Taille caractéristique des éléments du maillage en mètres
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", lc)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", lc)
# Génération du maillage volumique (3D)
gmsh.model.mesh.generate(3)
# Sauvegarde du maillage au format .msh
gmsh.write("ellipse.msh")
# Finalisation de GMSH
gmsh.finalize()
# Lecture du fichier .msh avec meshio
msh = meshio.read("ellipse.msh")
# Conversion et sauvegarde du maillage au format XDMF (compatible avec FEniCSx)
meshio.write("ellipse.xdmf", meshio.Mesh(points=msh.points, cells={"tetra": msh.cells_dict.get("tetra", [])}))
# Lecture du maillage XDMF avec FEniCSx
with XDMFFile(MPI.COMM_WORLD, "ellipse.xdmf", "r") as xdmf_file:
domain = xdmf_file.read_mesh(name="Grid")
domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)
# Affichage des informations sur le maillage
print("Maillage créé et importé avec succès.")
print(f"Nombre de sommets : {domain.topology.index_map(domain.topology.dim).size_global}")
print(f"Nombre d'éléments : {domain.topology.index_map(domain.topology.dim-1).size_global}")
Visualisation de la géométrie#
# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
interactive_panel = pn.pane.VTK(p.ren_win, width=700, height=400)
interactive_panel
from myst_nb import glue
glue("img", interactive_panel,display=False)
Configuration du problème#
import ipywidgets as widgets
from IPython.display import display
# Exemple de slider pour changer le rayon de la poutre
rayon = widgets.FloatSlider(value=0.05, min=0.01, max=0.1, step=0.01, description="Rayon")
display(rayon)
# Code Fenics ou autre qui dépend de "rayon.value"
import ipywidgets as widgets
from IPython.display import display
# Exemple de slider pour changer le rayon de la poutre
rayon = widgets.FloatSlider(value=0.05, min=0.01, max=0.1, step=0.01, description="Rayon")
display(rayon)
# Code Fenics ou autre qui dépend de "rayon.value"
Définition de l’espace fonctionnel#
# Exemple de dropdown pour choisir un matériau
materiau = widgets.Dropdown(
options=['Acier', 'Aluminium', 'Plastique'],
value='Acier',
description='Matériau:',
)
# Fonction pour mettre à jour le slider en fonction du matériau choisi
def update_slider(change):
if change['new'] == 'Acier':
rayon.value = 0.05 # Valeur par défaut pour l'acier
elif change['new'] == 'Aluminium':
rayon.value = 0.04 # Valeur par défaut pour l'aluminium
else:
rayon.value = 0.03 # Valeur par défaut pour le plastique
materiau.observe(update_slider, names='value')
imk=materiau.observe(update_slider, names='value')
display(materiau)
# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
interactive_panel = pn.pane.VTK(imk.ren_win, width=700, height=400)
from myst_nb import glue
glue("img", interactive_panel)
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Cell In[101], line 3
1 # Utiliser Panel pour créer un rendu interactif
2 pn.extension('vtk')
----> 3 interactive_panel = pn.pane.VTK(imk.ren_win, width=700, height=400)
4 from myst_nb import glue
5 glue("img", interactive_panel)
AttributeError: 'NoneType' object has no attribute 'ren_win'
Définition des frontière du domaine#
Show code cell content
# Définition des fonctions pour identifier les différentes faces du cylindre
def down_face(x):
"""Identifie la face inférieure du cylindre (z = 0)"""
return np.isclose(x[2], 0)
def top_face(x):
"""Identifie la face supérieure du cylindre (z = h)"""
return np.isclose(x[2], h)
def lateral_face(x):
"""
Identifie la surface latérale du cylindre elliptique
Utilise l'équation de l'ellipse : x^2/a^2 + y^2/b^2 = 1
"""
tolerance = 1e-5 # Tolérance pour la comparaison numérique
return (np.isclose((x[0]**2 / a**2 + x[1]**2 / b**2), 1.0, atol=tolerance)) & (0 <= x[2]) & (x[2] <= h)
# Dimension des facettes (2D pour un domaine 3D)
fdim = domain.topology.dim - 1
# Localisation des entités (facettes) correspondant à chaque face
Sigma_l = mesh.locate_entities(domain, fdim, lateral_face)
Sigma_0 = mesh.locate_entities(domain, fdim, down_face)
Sigma_h = mesh.locate_entities(domain, fdim, top_face)
# Préparation des marqueurs de facettes pour les conditions aux limites
marked_facets = np.hstack([Sigma_h]) # Utilisation uniquement de la face supérieure
marked_values = np.hstack([np.full_like(Sigma_h, 1)]) # Marqueur 1 pour la face supérieure
# Tri et création des tags de maillage
sorted_facets = np.argsort(marked_facets)
facet_tag = dolfinx.mesh.meshtags(domain, domain.topology.dim - 1,
marked_facets[sorted_facets],
marked_values[sorted_facets])
# Définition de la mesure pour les facettes marquées
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag)
Visualisation des frontière du domaine#
Show code cell outputs
# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
interactive_panel = pn.pane.VTK(pl.ren_win, width=600, height=400)
from myst_nb import glue
glue("img", interactive_panel)
Conditions aux limites#
# Définition des conditions aux limites de Dirichlet (déplacement nul)
u_D = np.array([0, 0, 0], dtype=default_scalar_type)
# Application des conditions aux limites sur chaque face
bc1 = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, Sigma_0), V) # Face inférieure
bc2 = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, Sigma_h), V) # Face supérieure
bc3 = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, Sigma_l), V) # Surface latérale
Propriétés matériau#
# Définition des constantes élastiques du matériau
mu = 80e9 # Module de cisaillement (G) en Pa
lambda_ = 120e9 # Premier paramètre de Lamé (λ) en Pa
Définition des chargements#
α = 20.0 m^-1
C = 1005309649.1487346 kg.m^2/s
Résolution#
Show code cell content
# Définition des tenseurs de déformation et de contrainte
def epsilon(u):
"""Tenseur de déformation linéarisé"""
return ufl.sym(ufl.grad(u))
def sigma(u):
"""Tenseur des contraintes (loi de Hooke)"""
return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)
# Définition des fonctions d'essai et de test
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
# Définition de la force volumique (nulle dans ce cas)
f = fem.Constant(domain, default_scalar_type((0, 0, 0)))
# Formulation variationnelle du problème
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx # Forme bilinéaire
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds(1) # Forme linéaire
# Définition et résolution du problème linéaire
problem = LinearProblem(a, L, bcs=[bc1],
petsc_options={"ksp_type": "cg", "pc_type": "jacobi"})
uh = problem.solve()
Visualisation des résultats#
Visualisation des déplacements#
# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
interactive_panel = pn.pane.VTK(p.ren_win, width=600, height=400)
from myst_nb import glue
glue("img", interactive_panel)
Visualisation des rotations principales#
Visualisation des déplacements amplifiés#
Visualisation des déplacements normalisés et selon les axes principaux#
Visualisation des composantes du tenseur des déformations#
Visualisation des composantes du tenseur des contraintes#
Calcul des contraintes de von Mises#
s = sigma(uh) - 1. / 3 * ufl.tr(sigma(uh)) * ufl.Identity(len(uh))
von_Mises = ufl.sqrt(3. / 2 * ufl.inner(s, s))
V_von_mises = fem.functionspace(domain, ("DG", 0))
stress_expr = fem.Expression(von_Mises, V_von_mises.element.interpolation_points())
stresses = fem.Function(V_von_mises)
stresses.interpolate(stress_expr)
Visualisation des contraintes de von Mises#
# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
interactive_panel = pn.pane.VTK(p.ren_win, width=600, height=400)
from myst_nb import glue
glue("img", interactive_panel)
Analyses des résultats#
Calcul des valeurs maximales#
Show code cell content
# Extraction de la valeur maximale de la rotation RZ
max_rotation_z = np.max(np.abs(rotation_z.x.array))
# Extraction de la valeur maximale de la norme de rotation
max_rotation_norm = np.max(rotation_norm.x.array)
# Conversion en degrés pour une interprétation plus intuitive
max_rotation_z_degrees = np.degrees(max_rotation_z)
# Calcul de l'angle de torsion théorique
theoretical_twist = ang * h # ang est l'angle de torsion par unité de longueur, h est la hauteur du cylindre
# Comparaison pour RZ
relative_error_Z = abs(max_rotation_z - theoretical_twist) / theoretical_twist * 100
# Comparaison pour norme R
relative_error_N = (max_rotation_norm - theta_theoretical) / theta_theoretical * 100
Moment d'inertie polaire J : 6.283185e-04 m^4
Angle de torsion théorique θ : 20.000000 radians
Angle de torsion théorique θ : 1145.92 degrés
Comparaison avec la simulation numérique en RZ:
Rotation Z maximale simulée : 19.888688 radians
Rotation Z maximale simulée : 1139.54 degrés
Erreur relative : 0.56%
Comparaison avec la simulation numérique en norme:
Norme de rotation maximale : 19.901232 radians
Norme de rotation maximale : 1140.26 degrés
Erreur relative : -0.49%
Valeurs maximales des composantes de rotation :
RX max : 124.69 degrés
RX max : 2.176298 radians
RY max : 122.15 degrés
RY max : 2.131950 radians